Draft version February 5, 2008 

Preprint typeset using I^Tf^X style emulateapj v. 6/22/04 



TWO EVOLUTIONAL PATHS OF AN AXISYMMETRIC GRAVITATIONAL 
INSTABILITY IN THE DUST LAYER OF A PROTOPLANETARY DISK 

FUMIHARU YAMOTO 1 
Division of Theoretical Astronomy, National Astronomical Observatory of Japan, 
Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan 

AND 

MlNORU Sekiya 
Department of Earth and Planetary Sciences, Faculty of Sciences, 
Kyushu University, Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan 
Draft version February 5, 2008 

ABSTRACT 

Nonlinear numerical simulations are performed to investigate the density evolution in the dust layer 
of a protoplanetary disk due to the gravitational instability and dust settling toward the midplane. 
We assume the region where the radial pressure gradient at equilibrium is negligible so that the shear- 
induced instability is avoided, and also restrict to an axisymmetric perturbation as a first step of 
nonlinear numerical simulations of the gravitational instability. We find that there are two different 
evolutional paths of the gravitational instability depending on the nondimensional gas friction time, 
which is defined as the product of the gas friction time and the Keplerian angular velocity. If the 
nondimensional gas friction time is equal to 0.01, the gravitational instability grows faster than dust 
settling. On the other hand, if the nondimensional gas friction time is equal to 0.1, dust aggregates 
settle sufficiently before the gravitational instability grows. In the latter case, an approximate ana- 
lytical calculation reveals that dust settling is faster than the growth of the gravitational instability 
regardless of the dust density at the midplane. Thus, the dust layer becomes extremely thin and may 
reach a few tenth of the material density of the dust before the gravitational instability grows. 
Subject headings: planetary systems: protoplanetary disks — solar system: formation — hydrody- 
namics — instabilities 



1. INTRODUCTION 

One of unsolved issues in the planetary system forma- 
tion is how kilometer-sized planetesimals form. There 
are two controversial models of the planetesimal for- 
mation: one is the growth through mutual stick- 
ing of dust aggregates due to nongravitational forces 
(iWeidenschilling 'fc Cuzzil 11993 iBIum fc WufTrl 2000: 
Sirono 2004), and another is the fragmentation of the 
dust layer due to the gravitational instability (G I) 
(|Safronovin963 lOoldreich fc Wa7dllT973 ISekivalll983h . 

Micron-sized dust particles in a protoplanetary disk 
grow due to the thermal motion in the initia l stage of the 
dust evolution (jHavashi fc Nakagawalll975t IB him et alJ 
2000). Unless the protoplanetary disk is globally turbu- 
lent, dust aggregates sweep up dust particles during the 
settling toward the midplane and the radial drift toward 
the central star, a nd gro w to centimeter-sized aggregates 
( Weidenschill ina 119801: iNakagawa. Sekiva. fc Havashil 
1986: iNomura fc Nakagawal \200ft) . When the dust 
density, which is defined as the total mass of the 
dust per unit volume of the disk, at the midplane 
exceeds a few tenth of the gas density, a velocity 
shear in the vertical direction grows. The shear- 
induced instability will occur due to the velocity 
shear and prevents dust settling ( Weidenschilling 1980; 
iCuzzi, Dobrovolskis, <fe ChampnevJ |199fl). Thus, the 
dust density at the midplane is difficult to approach the 
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critical density of the GI, which is several hundreds of 
the gas density. 

A few possibilities to exceed the critical density of the 
GI have been considered: (1) If dust aggregates grow 
to the radii of the order of 10 m, the coupling of the 
dust and gas becomes weak and the dust density at the 
midplane exceeds the critical density of the GI (Cuzzi e t 
al. 1993: iDobrovolskis. Dalles-Mar iani. fc Cuzzil 11999^. 
However, Weidenschilling (1995) has noted that the GI 
would be prevented by radial random velocities of the 
dust. (2) If the dust-to-gas surface density ratio in- 
creases by an order of magnitude through the dust 
migration or gas dis sipation, the shear-induced insta- 
bility is suppressed (ISekivl Il99?i |Youdin fc Shii|[2002: 
Youdin & Chiang] I2004T) . (3) If the pressure has a local 
maximum value in a protoplanetary disk, the radial pres- 
sure gradient which causes the shear-induced in stability 
would be negligible ( Haghighi pour fc Bossll2003i) . In this 
letter, we focus on the last case and other possibilities 
will be considered in future works. 

We have investigated the density evolution in the dust 
layer due to the GI and dust settling. In §2, a model of 
numerical simulations is described. In §3, results of nu- 
merical simulations and an approximate analytical calcu- 
lation are shown. Finally, these results are summarized 
in §4. 

2. MODEL 

We use the local Cartesian coordinate system (x, 
z) rotating with at a representative distance R from 
the central star, where Qk is the circular Keplerian angu- 
lar velocity, instead of the inertial system (r, </>, z) with 
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(a) t = 0.0 (b) 7=3.0 
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Fig. 1. — The density evolution by the GI and dust settling for if = 0.01. Figures show the snapshots of the density at (a) t = 0.0 (i.e., 
the initial condition), (b) t = 3.0, (c) t = 5.5, and (d) t = 6.3. The length scales in the x and z directions are normalized by the scale 
height of the dust layer h&. The gas is incompressible and has a constant density. I n this simulation, we assume p c /pg = 340. This value 
corresponds to the gas density at 1 AU for the minimum mass solar nebula (Havashi 1981), for which we have = 86 km. The numerical 
simulation is performed for the region between the midplane (z = 0) and the upper boundary (z = 4.0/id), and lower halves of figures are 
drawn by using mirror symmetry with respect to the midplane. 



its origin at the central star. We use approximations in 
which only linear terms of (x, y, z) are left and higher 
order terms are neglected: r ~ R + x and y ~ R<j). In this 
letter, the dust layer is assumed to be axisymmetric with 
respect to the rotation axis, i.e., any physical quantity is 
independent of y. We assume the region where the radial 
pressure gradient at equilibrium is negligible so that the 
shear-induced instability is avoided. 

The equations of continuity and motion of the dust and 
gas, and the Poisson equation for the self-gravitational 
potential are solved numerically. The tidal, the Coriolis 
and the self-gravity forces on the dust and gas are ap- 
plied. Additionally, the gas drag force on the dust and 
its reaction and the pressure gradient forces on the gas 
are applied. The gas drag and its reaction forces de- 
pend on the gas friction time tf, which is defined as the 
time during which the relative velocity between a dust 
aggregate and gas becomes 1/e. The dust pressure gra- 
dient force is assumed to be neglected because random 
velocities of dust aggregates are sufficiently suppressed 
by the gas drag force before the GI grows if U <^ ^k*- 
All dust aggregates are assumed to have the same gas 
friction time which is constant for time. The gas is as- 
sumed to be incompressible. This assumpti on is good in 
the dust layer as explained in iSekival (|1983t ). 

The Gaussian dust density distribution in the direc- 
tion vertical to the midplane is used as the unperturbed 
density. The initial value of the density at the midplane 
is assumed to be equal to the critical density of the GI p c 



(jYamoto fc Sekival 120041) . Unperturbed velocities of the 
dust and gas are set to zero in the x direction and the 
Keplerian velocity in the y direction because we assume 
that the radial pressure gradient is negligible. In the z 
direction, the unperturbed gas velocity is set to zero and 
the unperturbed dust velocity is given by the terminal 
settling velocity: 

w d = -tf^Q^z + AttG J (p g + Pd)dz}, (1) 

where p g and pa is the density of the gas and dust, re- 
spectively. 

As an initial condition, we give only the dust density 
perturbation which has the criti cal wavelength of the GI 
obtained by the linear analysis (jYamoto fc Sekival l2004) 
and the maximum amplitude of a hundredth of the un- 
perturbed dust density. The perturbed velocities of the 
dust and gas are set to be zero. In numerical simulations, 
boundaries in the x direction are assumed to have a pe- 
riod of the critical wavelength. We assume that the gas 
density and the dust surface density are equal to those 
of the minimum mass solar nebula mod el (Havashi 1981; 
lHayashi, Nakazawa, <fe Nakae?awa|l98^ at 1 AU, but re- 
sults are not sensitive to these assumptions and the only 
important parameter is the nondimensionalized gas fric- 
tion time tf (hereafter we denote nondimensional time 
and density by diacritic tilde, which are normalized by 
Q^ 1 and Q|-/(47rG) respectively). 

3. RESULTS 



Gravitational instability in the dust layer 
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The density evolution for tf = 0.01 is shown in Fig- 
ure H I n the early stage of the numerical simulation 
(see Figs. and^), the growth rate of the GI is so 
small that the density around the midplane mainly in- 
creases by dust settling. Next, the GI grows faster than 
dust settling and the dust layer in the calculated region 
changes into a ring circling the central star (see Fig. QJ:). 
Afterwards, the density has a maximum value on the left 
and right in the ring and concave density distributions 
arise in the vicinity of a maximum density (see Fig. Qi). 
Subsequent evolution of the density is significant, but we 
should perform three dimensional numerical simulations 
because nonaxisymmetric modes may grow in this stage. 

The snapshot of the density at t = 0.9 for tf = 0.1 
is shown in Figure The settling velocity of dust ag- 
gregates (eq. [Q) for tf = 0.1 is ten times as large as 
that for tf = 0.01 due to the difference of the gas friction 
time. The numerical simulation shows that dust settling 
is faster than the growth of the GI as long as the dust 
density at the midplane is smaller than about twenty 
times of the critical density of the GI. 
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Fig. 2. — The snapshot of the density at t = 0.95 for tf = 0.1. 
The length scales, the gas density, the upper boundary, and the 
simulated region are the same as Figure^ Note that the maximum 
value of the density is different from Figure [I] 



For t{ = 0.1, dust settling is faster than the growth of 
the GI within the time interval of the numerical simu- 
lation. However, the numerical simulation is performed 
for the density only up to about twenty times the critical 
density. Here, we show that dust settling is faster than 
the growth of the GI regardless of the dust density at the 
midplane. The dust settling time Nettie during which the 
dust density increases by a factor of e is given by 

(1 + p g )tf [ V Pa + 1 + Pg J J 

using a rough approximation described in Appendix A. 
If pa ^> 1 + p g : the vertical component of the self-gravity 
of the dust in the dust layer is much larger than the 
sum of vertical components of the central star's gravity 
and the self-gravity of the gas, the dust settling time is 
approximately equal to [1 — (1 / e)] / (pdU) . On the other 

hand, if p d < 1 + p g , Settle ~ VK 1 + Pg)tf], and Settle 

results in 1/tf for the minimum mass solar nebula model 
because p g <C 1. Figure [3| shows the dust settling time 
with t f = 0.01 and 0.1. 



Fig. 3. — Comparison between the dust settling time toward 
the midplane with tf = 0.01 (upper solid line) and 0.1 (lower solid 
line) and the growth time of the GI with the Gaussian dust density 
distribution (bold line). The gas density p g is shown by the dotted 
line. The critical density of the GI p c is shown by the dashed line. 



The growth time of the GI is given by the reciprocal 
of its growth rate. The growth rate is obtained using 
the m ethod written in §3 and §4 of lYamoto fc Sekival 
(2004) by assuming that the dust density distribution 
in the direction vertical to the midplane is kept to be 
Gaussian regardless of the value of the dust density at the 
midplane. Although the dust density distribution will 
diverge from the Gaussian as dust aggregates settle, this 
assumption is fairly good for our rough estimate. The 
growth time of the GI with the Gaussian dust density 
distribution is shown in Figure 01 for comparison with 
the dust settling time. After the dust density exceeds 
the critical density, the growth time of the GI decreases 
rapidly, and is proportional to the dust density to the 
power —1/2 for pd ^> p c . On the other hand, the dust 
settling time is the constant value for pd < p g , and is 
proportional to the dust density to the power —1 for 
Pd ^> Pg- Figure 01 shows that the GI grows faster than 
dust settling up to thirty times of the critical density for 
tf = 0.01; on the other hand, dust settling is always faster 
than the growth of the GI for tf = 0.1. The results of 
the approximate analytical calculation agree with those 
of the numerical simulations. 

4. SUMMARY 

We have investigated the density evolution in the dust 
layer due to the GI and dust settling by numerical sim- 
ulations and the analytical calculation. If tf = 0.01, the 
numerical simulation shows that the GI grows faster than 
dust settling. This path of the GI is con sidered to be 
analogous with the case of iSekival ((1983) applying the 
one-fluid model which the dust and gas have the same 
velocity due to the strong coupling by the gas drag force, 
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although our simulations take into account the velocity 
difference of the dust and gas. On the other hand, if 
ti — 0.1, the numerical simulation shows that dust set- 
tling is faster than the growth of the GI as long as the 
dust density at the midplane is smaller than about twenty 
times of the critical density of the GI. In this case, the 
approximate analytical calculation reveals that dust set- 
tling is faster than the growth of the GI regardless of 
the dust density at the midplane, so the dust layer will 
become extremely thin. Subsequently, the GI will grow 
in the extreme thin dust layer with much smaller wave- 
lengths than the critical wavelength, and then modes 
with larger wavelengths may begin to grow slowly. This 
p ath of the GI is considered to be analogous with the case 
of lGoldreich fc \fy ard (1973). Therefore, our results indi- 
cate that the difference of the nondimensional gas friction 
time leads to two evolutional paths of the axisymmetric 
GI. 



In order to confirm whether the planetesimal formation 
by the GI is possible or not, nonaxisymmetric calcula- 
tions are essential. We plan to perform the numerical 
simulations and investigate three dimensional evolution 
of the GI in subsequent papers. 



Numerical simulations have been performed on 
VPP5000 at the Astronomical Data Center of the Na- 
tional Astronomical Observatory of Japan and on eS- 
erver p5 595 at Computing and Communications Center 
of Kyushu University. F. Y. is supported by the Re- 
search Fellowship of Japan Society for the Promotion of 
Science (JSPS) for Young Scientists. M. S. is supported 
by Grants-in-Aid of the Japanese Ministry of Education, 
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APPENDIX 

DERIVATION OF EQUATION (2) 

Here, we show how to derive equation (0). Although the dust density is not constant for z, we here assume dpa/dz = 
and d(pdUd) / dx = for a rough estimate of the growth time of the dust density. Then the continuity equation of the 
dust is written 



dpd 
dt 



Pd- 



dwd 
dz 



0. 



Substituting equation into equation (|A1|) and using nondimensional variables, we have 

dpd 



dt 

Integrating equation (|A2|) with respect to f, we have 

1 



pdU(pd + 1 + p g ). 



t 



-^ln 



(l+p g )tf VPd + 1 + 



Pd 



(Al) 



(A2) 



(A3) 



where C is an arbitrary constant of integration. From equation (|A3|) . the dust settling time during which the dust 
density increases by a factor of e results in equation (J2J). 
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